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Abstract 


We  consider  the  use  of  fast  direct  methods  as  preeonditioners  for  iterative  methods  for  computing  the 
numerical  solution  of  non- self-adjoint  elliptic  boundary  value  problems.  We  derive  bounds  on 
convergence  rates  that  are  independent  of  discretization  mesh  size.  For  two-dimensional  problems  on 
rectangular  domains,  discretized  on  an  nin  grid,  these  bounds  lead  to  asymptotic  operation  counts  of 
0(n2log  n  log  tl)  to  achieve  relative  error  <  and  0(n2(log  n)2)  to  reach  truncation  error. 
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1.  Introduction 


The  numerical  solution  of  elliptic  boundary  value  problems  by  finite  difference  methods  requires  the 
solution  of  systems  of  linear  equations  of  the  form 

Au  =  f,  (1) 

where  A  is  a  large  sparse  matrix.  Fast  direct  methods  are  efficient  techniques  for  solving  (1)  when  the 
elliptic  equation  is  posed  on  a  rectangular  domain  and  the  differential  operator  is  separable.  In  this 
paper,  we  show  that  nonseparable,  non-self-adjoint  elliptic  problems  can  be  solved  efficiently  using  fast 
direct  methods  as  preconditioners  for  iterative  methods  for  nonsymmetric  linear  systems. 

As  a  prototype,  consider  the  two-dimensional  Dirichlet  problem 


Au  =  /,  u  G  0,  (2) 

u  —  g,  u  €  3/7, 


where  fJ  is  a  rectangular  region  in  R2, 


Au  m  -  (aux)x  -  (buy)y  +  c«x  +  (cu)x  +  duy  +  (du)y  +  eu,  (3) 

and  a(x,y),  b(x,y),  c(x,y),  d(x,y),  e(x,y),  and  /(x,y)  are  smooth  functions  defined  on  fi,  with  a,  b  >  0, 
e  >  0.  Discretizing  (2)  by  finite  difference  techniques  on  an  nxn  grid  leads  to  a  sparse  system  of  linear 
equations  of  the  form  (1)  where  A  is  of  order  N  =  n2.  If  c  and  d  are  zero,  then  A  is  self-adjoint; 
otherwise,  A  is  non-self-adjoint  and,  in  general,  A  is  nonsymmetric.  If 


a  —  a(x),  b  —  b(y),  c  =  c(x),  d  =  d(y),  e  —  e^x)  +  e2(y), 


(4) 


then  A  is  separable.  In  the  self-adjoint  separable  case,  (1)  can  be  solved  by  a  variety  of  fast  direct 
methods,  such  as  the  cyclic  reduction  and  Fourier  methods  (surveyed  in  [7,  18])  and  the  generalized 
marching  algorithm  [2,  3].  In  the  non-self-adjoint  case,  the  cyclic  reduction  method  can  still  be  used  (17]. 
All  of  these  methods  require  0(n2log  n)  arithmetic  operations  (i.e.  floating  point  multiplications  and 
divisions).1 

If  A  is  nonseparable,  then  fast  direct  methods  can  be  used  to  solve  (1)  iteratively.  For  self-adjoint 
problems,  Widlund  [23]  proposed  the  stationary  method 


ui+l  “  ui  +  rQ',(f*Aui). 


’in  the  self-adjoint  case,  this  count  can  be  reduced  to  OJn^ogflng  n)).  See  (18]. 


/ 


(5) 
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where  Q  is  the  discretization  of  a  self-adjoint  separable  approximation  of  A.  Concus  and  Golub  [5]  and 
Bank  [2]  examined  accelerating  (5)  by  the  Chebyshev  and  conjugate  gradient  methods  respectively; 
equivalently,  Q  is  a  preconditioning  matrix  for  these  iterative  methods.  Convergence  of  all  these 
methods  is  independent  of  the  mesh  size  used  in  the  discretization,  so  that  a  relative  error  of  t  is 
achieved  in  0(!og  e1)  iterations.  Since  the  dominant  expense  of  each  iteration  is  a  fast  direct  solve,  the 
asymptotic  operation  count  is  0(n2log  n  log  c*1). 

Concus  and  Golub  (4]  and  Widlund  [22]  extended  this  analysis  to  some  particular  non-self- adjoint 
problems  using  the  generalized  conjugate  gradient  method,  which  depends  on  the  symmetric  part  of  A  as 
a  preconditioner.  Let 

Mu  m  -  (atijj  -  (buy)y  +  eu,  (6) 

Bu  2S  cux  +  (cu)K  +  duy  +  (du)y, 

so  that  A  “  M  +  R.  The  symmetric  part  of  A  corresponds  to  the  discretization  of  M,  the  sum  of  the 
second  and  zero  order  terms  of  A.  The  convergence  of  the  generalized  conjugate  gradient  method 
depends  essentially  on  the  spectral  radius  of  the  discrete  analogue  of  the  compact  operator  MlR,  so 
that  convergence  is  independent  of  mesh  size.  If  M  is  separable,  i.e.  a,  b,  and  e  are  as  in  (4),  then  fast 
direct  methods  can  be  used  for  the  preconditioning,  so  that  the  asymptotic  operation  count  is  again 
0(n2log  n  log  c'1). 

In  this  paper  we  show  that  these  asymptotic  operation  counts  can  be  achieved  even  if  the  symmetric 
part  of  A  does  not  come  from  a  separable  operator.  Other  iterative  methods  for  no nsym metric  linear 
systems,  such  as  the  conjugate  gradient  method  applied  to  the  normal  equations  [13]  and 
Orthomin(k)  [8,  9],  allow  more  general  choices  than  the  symmetric  part  for  preconditioning  matrices. 
Using  an  analysis  of  the  finite  difference  discretization,  we  show  that  with  symmetric  positive-definite 
preconditioning  matrices  derived  from  other  self-adjoint,  separable  operators  Q  that  approximate  A,  the 
asymptotic  convergence  rates  of  these  iterative  methods  is  independent  of  mesh  size.  Although  this 
analysis  holds  only  for  symmetric  preconditioning  operators,  we  provide  numerical  examples  that  suggest 
that  nonsymmetric  preconditioners  derived  from  non-self-adjoint  separable  approximations  of  A  lead  to 
the  same  asymptotic  convergence  properties.  In  Section  2  we  review  the  basic  properties  and 
convergence  bounds  of  the  conjugate  gradient  method  applied  to  the  normal  equations  and  Orthomin(k). 
In  Section  3,  we  present  the  convergence  analysis  for  general  self-adjoint,  separable  preconditioning 
operators,  and  in  Section  4,  we  demonstrate  the  performance  of  these  techniques  with  both  self-adjoint 


and  non-self-adjoint  preconditioning  operators. 


2.  Convergence  Bounds  for  Iterative  Methods 

Given  a  nonsymmetric  linear  system  of  the  form  (1),  let  A  =  M  +  R,  where  M  —  (A+AT)/2  is  the 
symmetric  part  of  A  and  R  ==  (A-A^)/2  is  the  skew-symmetric  part  of  A.  In  this  section,  we  present 
upper  bounds  on  the  convergence  rates  of  iterative  methods  for  nonsymmetric  linear  systems  in  terms  of 
eigenvalues  of  M  and  R.  We  consider  two  iterative  methods:  the  conjugate  gradient  method  applied  to 
the  normal  equations  and,  as  representative  of  a  recently  developed  collection  of  methods  that  avoid  the 
use  of  the  normal  equations,  Orthomin(k). 

We  first  establish  some  conventions  of  notation.  For  any  square  matrix  B,  let  ic(B)  as  ||B||2||B',||2 
denote  the  condition  number  of  B,  let  X  .  (B)  and  (B)  denote  the  eigenvalues  of  B  with  smallest  and 
largest  modulus,  respectively,  and  let  />(B)  denote  the  spectral  radius  |XmM(B)|.  For  symmetric  positive- 
definite  B,  let  j|v||B  denote  the  B-norm  (v.Bv)1/2. 

A  classical  method  for  solving  nonsymmetric  linear  systems  of  the  form  (1)  is  to  apply  the  conjugate 
gradient  method  (CG)  [13]  to  the  normal  equations 

ATAu  =  ATf. 

That  is,  the  nonsymmetric  system  is  embedded  into  a  system  with  a  symmetric  positive-definite 
coefficient  matrix,  to  which  the  conjugate  gradient  method  is  applicable.  We  denote  this  method  by 
CGN.2  The  convergence  properties  of  CGN,  derived  from  the  standard  error  analysis  of  CG  [0],  arc  well 
understood.  We  summarize  the  results  we  need  as  follows.  Let  {u;}  denote  the  sequence  of  iterates 
generated  by  CGN  starting  from  an  initial  guess  uQ,  and  let  {rj  denote  the  residuals  {f-Auj}.  The 
iterate  u;  is  the  point  in 

uQ  +  span  {ATr0,(ATA)ATr0,...,(ATA)l  IATr0} 
whose  residual  norm  [|rj||2  is  minimum.  As  a  result,  the  residuals  satisfy  [6] 

»*  2  2  Mr  PI 

The  following  result  relates  this  bound  to  the  extreme  eigenvalues  of  M  and  the  spectral  radius  of  R. 


2A  related  technique  is  to  solve  AATu'««f  by  CG,  with  U“ATu'. 
properties  as  CGN,  see  [9], 


This  method  has  essentially  the  same  convergence 


Theorem  1:  If  the  symmetric  part  M  of  A  is  positive-definite,  then 


*m;„<ATA)  >  X  .  <M)J, 


>W(ATA)  <  |X„„(M)  +  ^R)| 


Hence,  the  residuals  {r;}  generated  by  CGN  satisfy 

■'A  - 2 !' '  owfew  + 1  ]' 


max'  '  r'  min' 

Proof:  Let  S  denote  the  unique  symmetric  positive-definite  square  root  of  M,  i.e.  S2  =  M. 
Then 

(v,ATAv)  —  (Av.Av)  -  (SlS+S^RJv.StS+S-'RJv) 

—  ((S+S^RJr^S+S-’RW. 

But  for  any  real  w, 

(*,Mw)  >  Xmin(MXw,w) 
and 

(w.Rw)  =  0, 
so  that 

(v,ATAv)  >  Xmin(M)((S+S-,R)v,(S+S-,R)v) 

-  Xmin(M)  [(Sv,Sv)  +  2(v,Rv)  +  (S^Rv.S^Rv)] 

=  Xmin(M)  [(v,Mv)  +  (S-,Rv,S-,Rv)) 


>  Xmin(M)  (v,v) 


Therefore, 

X  .  (AtA)  =  min  X  .  (M)2. 

min'  ’  T^0  (v,v)  ~  min'  ’ 

m 

For  the  upper  bound  on  XmM(A*A), 

Xm„(ATA)  -  ]|A|||  <  |(M|,  +  ||R|y2  -  [Xm„(M)  +  *R)|2, 
where  we  have  used  the  fact  that  ||R||2  =  p  since  R  is  skew-symmetric  and  hence  normal  [20). 


Finally,  note  that  the  fraction  in  (7)  satisfies 


The  bound  (7)  implies  that  CGN  is  convergent  for  arbitrary  nonsingolar  linear  systems,  and  as  we  will 
show  in  Theorem  5,  the  corollary  result  (10)  gives  rise  to  strong  asymptotic  bounds  under  suitable 
conditions.  However,  the  dependence  of  CGN  on  ATA  is  a  drawback.  If  A  itself  were  symmetric 
positive-definite,  then  CG  could  be  applied  directly  to  (1).  In  this  case,  the  bound  on  the  relative  error 
of  the  i’th  iterate  generated  by  CG  is  proportional  to 

MW- 

(see  [8]),  and  the  convergence  of  CG  would  be  much  faster  than  the  convergence  of  CGN.  Moreover,  to 
avoid  the  explicit  computation  of  A  A,  CGN  requires  two  matrix-vector  products  per  iteration,  one  by  A 
and  one  by  AT.  In  contrast,  CG  requires  only  one  matrix-vector  product  per  iteration  for  symmetric 
positive-definite  systems. 

In  recent  years,  many  conjugate  gradient-like  iterative  methods  for  nonsymmetric  systems  have  been 
developed  with  the  aim  of  avoiding  these  difficulties  of  CGN  (1,  8,  9, 16,  21,  25).  In  all  of  these 
methods,  the  i’th  iterate  u;  is  chosen  from  the  Kryiov  space  uQ  4-  K~,  where 

K;  »  span{r0,Ar0 . 

with  the  hope  that  the  convergence  depends  on  k(A)  rather  than  k(ATA).  Each  iteration  requires  one 
matrix-vector  product  of  the  form  Av.  In  this  paper,  we  will  take  the  method  known  as  Orthomin(k)  to 
be  representative  of  this  class  of  iterative  methods.  Orthomin(k)  chooses  for  u;  a  point  in  ufl  +  Kj  whose 
residual  norm  Ijrjflj  i»  minimum  in  a  (k+1  ^dimensional  subspace  of  K(  (see  [8]).  Although  its 
convergence  properties  are  not  as  well  understood  as  those  of  CG,  the  following  result  shows  that  when 
M  is  positive-definite,  Orthomin(k)  generates  a  sequence  of  approximate  solutions  that  converge  to  A’*f. 

Theorem  2t  The  residuals  (r.)  generated  by  Orthomin(k)  satisfy  [8,  9] 


K  JM>  +  f(R)2Ami„(M) 


The  work  per  step  of  the  two  methods  considered  in  this  section  is  [8,  9] 


CGN:  5N  operations  plus  two  matrix-vector  products,  Av  and  ATv; 

Orthomin(k):  (3k+4)N  operations  plus  one  matrix-vector  product,  Av. 


3.  Convergence  Analysis  for  Separable  Preconditioners 

In  this  section,  we  »;how  that  if  CGN  or  Orthomin(k)  is  combined  with  preconditioners  based  on 
certain  separable  operators  to  solve  discretized  elliptic  problems,  then  convergence  does  not  depend  on 
the  mesh  size  used  in  the  discretization. 


Let  A,  /  and  g  be  defined  as  in  (2),  let  fl  denote  the  unit  square  0<x,y<l,  and  let  A  »  M  +  R,  as 
in  (6).  Discretizing  (2)  by  centered  differences  on  an  nxn  grid  leads  to  a  system  of  linear  equations  of  the 
form  (1).  Let  h  =  l/(n+l).  In  terms  of  the  contributions  of  the  symmetric  and  skew-symmetric  parts, 
the  difference  equations  at  a  typical  mesh  point  (xj,yj)  have  the  form  (after  scaling  by  h2) 

(Au]"  —  [Mu]-  +  [Ru]-  —  h2f;j  , 


where 


[Mu]-  -  [ai+i/2,j+ai-!/2,j  +  bi,j+l/2+bi,j+l/2  +  h2eijluij 

*  ai+l/2,jui+I,j  '  a*-l/2,jUi-l j  '  bi,j+l/2Ui,j+l  '  bi,j-l/2ui,j-l> 

IHj  =  Kci+l,j+cij)h/2]ui+l,j  *  KCij+CH,i)h/2lUi-l,j  +  Kdi,j+l+dij)h/2lui,j+l 

*  Kdij+di,j-l)h/2lui,j-|- 

Hence,  M  corresponds  to  the  discretization  of  M,  and  R  corresponds  to  the  discretization  of  R. 


(12) 


(13) 


Let  Q  denote  a  self-adjoint  elliptic  operator  on  /?  which  is  spectrally  equivalent  to  M,  i.e. 


for  all  sufficiently  smooth  v  t  L2(/7)  satisfying  homogeneous  Diricblet  boundary  conditions,  where 

(»»,«;)  s*  /  wu> , 

J  n 

and  o,  and  are  positive  constants.  For  example,  Q  could  be  the  negative  Laplacian  4-  — J. 

1  2  Idx2  dy2J 

The  discretization  of  Q  gives  rise  to  a  symmetric  positive-definite  matrix,  Q,  that  is  spectrally  equivalent 

to  M:  a.  and  q2  in  (14)  can  be  chosen  so  that 


independent  of  the  mesh  site  used  in  generating  M  and  Q. 


Since  Q  is  symmetric  positive-definite,  it  admits  a  factorisation  Q  =*  LLT  (here  L  is  not  necessarily 
lower  triangular).  Consider  the  eymmttrieally  preconditioned  linear  system 

Au  —  [L**AL  T)  (LTu]  =  L  *f  =  f  (16) 

in  which  the  coefficient  matrix  A  has  symmetric  part  M  =  L**ML'T  and  skew-symmetric  part  R  = 
L'*RL'T.  The  extreme  eigenvalues  of  M  are  given  by 


X  ■  (M)  -  min  (y^MLfyjL 
mm'  •  (v,v) 


Xm  JM)  =  max  faJ^ML'Trl 
T^o  (v,v) 


But  for  any  vy^O,  let  w  =  L'^v,  so  that 


fv.L^ML^vl  _  (w.Mw) 
(v,v)  KQw) 


(17) 


Hence,  by  (15), 

>  «|.  Xmmx(M)  <a2, 

independent  of  mesh  site. 


(18) 


Observation  (18)  is  the  basis  of  the  convergence  results  for  the  self-adjoint  case  [5].  To  extend  the 
analysis  to  the  non-self- adjoint  case,  it  is  necessary  to  bound  p{R).  An  intuitive  approach  is  to  note  that 
R  is  similar  to  Q'*R,  which  is  the  discrete  analogue  of  the  continuous  operator  T  s  Q'lR.  For  Q  with 
sufficiently  smooth  coefficients,  Tis  a  compact  operator  whose  eigenvalues  are  therefore  bounded  with  0 
as  their  sole  accumulation  point  [22,  24];  one  would  expect  a  similar  bound  on  the  eigenvalues  of  the 
discrete  operator.  We  derive  such  a  bound  below.  For  simplicity,  we  present  the  analysis  for  the  one¬ 
dimensional  Dirichlet  problem  on  the  unit  interval  Cl  [0,1],  divided  into  n  uniformly  spaced  interior 
mesh  points,  with  h=l/(n+l).  The  analysis  for  the  two-dimensional  Dirichlet  problem  on  a  rectangular 
domain  is  identical,  using  the  two-dimensional  versions  of  the  operators  in  (21)  below. 


The  one-dimensional  analogue  of  the  first  order  operator  R  in  (6)  is 
Ru  3  cux  +  (ctt)x  , 

and  the  discretization  analogous  to  that  of  (13)  is 

IRu)i  "  ((ci+i+ci)h/2]ui+1  -  [(ci+cM)h/2]ui.1  -  |  MVfUn)  +  ci+iui+rci-iui-il> 


(19) 


1  <  i  <  n. 
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Without  loss  of  generality,  we  take  Q  to  be  the  one-dimensional  Laplacian,  Qu  s  - u (The  argument 
can  be  modified  easily  for  any  operator  spectrally  equivalent  to  this  choice  of  Q.)  The  discretization 
analogous  to  that  of  (12)  is 

[Qu];  m  -ui+,  +  2uj  -  Uj.,,  1  <  i  <  n.  (20) 

In  (19)  and  (20),  u;,  [RuJj  and  (Qu]j  are  defined  to  be  zero  for  i— 0  and  i=n+l. 


Let  the  usual  finite  difference  operators  be  given  by 

0<i<D.  |D+.|nt,  -  0. 

u.-u-  .  (21) 

(D  u).  =  1  <  i  <  n+1,  (D_u]0  =  0. 

Note  that  [D+u]0  and  [D  u]n+1  may  be  nonzero  even  if  u0  =  un+1  =  0. 


To  avoid  repeatedly  handling  the  0’th  and  (n+l)’st  indices  as  special  cases,  we  adopt  the  following 
convention  of  notation.  We  identify  Cn  with  the  proper  subspace  of  Cn+2  consisting  of  complex  (n+2)- 

vectors  v  =  (v^ . v„.vn+i)T  such  that  vo  =  vn+i  “  Thus,  R  and  Q  of  (19)  and  (20)  are  defined 

on  this  representation  of  Cn  without  reference  to  indices  0  and  n+1.  The  usual  Euclidian  inner  product 
and  norm  on  Cn  are  inherited  by  this  space: 

(v,w)=  EViwi(  ||v||2  e  (v,v)'/2.  (22) 

We  will  also  abuse  notation  slightly  and  define  (v,w)  and  ||v||2  as  in  (22)  for  arbitrary  vectors  in  Cn+2, 
noting  that  these  functions  do  not  define  an  inner  product  or  norm  on  all  of  Cn+2.  We  refer  to  the  true 
Euclidian  inner  product  on  Cn+2  as  the  “extended”  inner  product 

■  +I 

(u,v)  S3  E  VjWj. 

i^O  1  1 

The  following  result  summarizes  the  properties  of  the  finite  difference  operators  that  we  need.  Other 
results  of  this  type  can  be  found  in  [14]. 

Lemma  3: 

(i)  For  v  and  w  €  Cn+2,  (v,D+w)e  -(D  v,w)e  +  J;!*n+1*n+,  -  v^J; 

For  v  and  w  €  Cn, 


(»»)  KD+w)  =  (v,D+w)e  =  -(D  v,w)e  =»  -(D  v,w); 

(iii)  Qv  =  -h2D  D+v  =  -h2D+D  v; 

(iv)  (v,Qv)  >  h2j|D+v||2  ,  (v,Qv)  >  h2||D  v|j2  . 


j  .v  --'--“v 


c.  <: 


Proof:  For  (i), 


WV->.  =  |vZi±p  _  -  E,i#i 


i+l 


I  »♦!  1*1  - 

“Kf,vi-iwi  '  gviwil  +  E  K+i*„+i  -  vowol 

—  *  (D_v,w)e  +  ifv.Vi-VoiJ 

Assertion  (ii)  follows  immediately  from  (i)  and  the  rcro  boundary  conditions. 
For  the  first  equality  of  (iii), 

IMrlMi-i  i  fv,+i-vi  vrvi-i 


The  proof  of  the  second  equality  is  identical. 


For  the  first  assertion  of  (iv), 


(v.Qv)  =*  (v,Qv)e  by  tb 

—  -h2(v,D  D+v)e  by  (ii 

—  h2(D+v.D+v)e  by  th 

>  h2||D+v||2  . 

The  proof  of  the  second  assertion  of  (iv)  is  identical. 


by  the  first  equality  of  (ii) 
by  (iii) 

by  the  second  equality  of  (ii) 


t  R  =  L'iRL'T,  where  Q  =  LLT  as  above. 

Theorem  4:  There  exists  a  constant  /?  >  0,  independent  of  the  mesh  sire  h,  such  that 
dQ-'R)  -  ,(R)  <  f). 

Proof:  Since  R  is  skew-symmetric  and  therefore  normal,  its  spectral  radius  is  given  by  the 
maximum  value  of  the  Rayleigh  quotient: 


v.L-’RL^v 


»ecr»^o  ii*,* 


Writing  Rw  in  terms  of  the  finite  difference  operators  of  (21), 
Rw  -  |i[c  (D++D>  +  (D++DJ(c  w)l  , 
where  c-w  denotes  the  vector  with  entries  {CjW.}^^.  Therefore 


(w,Rw)  =  jj-  [  (w,c  (D++D  )w)  +  (w,(D++D  )(c  w))  ] 

=  [  (c  w,(D++D  )w)  -  ((D++D>,(c'»))], 

by  Lemma  3,  (ii).  Applying  Cauchy-Schwari  and  letting  C  =  max  |c(x)|, 

i  en 

|(w,Rw)|  <  h2  ||c  w||2  ||(D++D_)w||2  <  C  h2  ||w||2(||D+w||2  +  ||D  w||2) 
<  2Ch  ||w||2  (w.Qw)1/2, 
by  Lemma  3,  (iv). 


^  „  hHI2  ^  a 

<  2C  —\j  <  0, 
(w,Qw)  '2 


with  0  independent  of  h,  since 

>  Xmi0(Q)  -  0<hJ)  (101. 


C.f.  (15]  for  an  analysis  of  the  convergence  of  the  eigenvalues  of  Q_1R  to  those  of  Q  !R. 
We  now  return  to  the  consideration  of  two-dimensional  problems. 

Note  that  applying  either  CGN  or  Orthomin(k)  formally  to  (16)  results  in  the  residuals 


Q.E.D. 


r  •  =  f  -  Au-  =  L'*r , 


so  that  the  residual  norm  associated  with  (16)  is  (IL^rJIg  =  IItjIIq-i-  Hence,  combining  (15),  (18)  and  (23) 
with  the  results  of  Theorems  1  and  2: 

Theorem  6:  If  Q  is  a  symmetric  matrix  that  satisfies  (15),  (18)  and  (23),  then  the  residuals 
generated  by  CGN  with  symmetric  preconditioning  by  Q  satisfy 


^  1’ 


and  the  residuals  generated  by  Orthomin(k)  with  symmetric  preconditioning  by  Q  satisfy 


Mq-i  <  1 


a,a2+^2 


i/2  Mo-.  • 


Hence,  for  discretized  elliptic  problems,  the  number  of  iterations  needed  to  make  ||r||n-i/llrnlln-i 


11 


<  <  is  independent  of  mesh  size. 

Proofs  For  the  first  inequality,  by  (18)  and  (23) 

(Xmu(M)+,KR))/Xmi[i(M)  <  («2+W«,. 
Hence,  applying  Theorem  1, 


M<*'  s  Ml  -  |„i+/)|/a|  + 


2 


°2*°1+^ 

02"*"ai"^ 


lfr0!lQ  i  • 


The  second  inequality  is  proved  in  a  similar  manner  using  Theorem  2. 


Q.E.D. 


The  extra  work  per  step  required  for  preconditioning  comes  in  the  matrix  vector  products  L"*AL*Tv 
and  (L'*AL'TjTv.  Preconditioned  CGN  and  Orthomin(k)  can  be  implemented  so  that  the  only  references 
to  Q  have  the  form  of  solving  Qw  =  v,  so  that  no  factorization  of  Q  is  required  (see  [9]).  The  solve 
requires  0(n2log  n)  operations,  but  the  remaining  work  per  step  of  both  iterative  methods  is  O(N)  = 
0(n2).  Hence,  asymptotically,  the  preconditioning  solves  are  the  dominant  cost  per  step. 

Corollary  6:  The  number  of  arithmetic  operations  required  by  either  CGN  or  Orthomin(k) 
with  preconditioning  by  Q  to  solve  (1)  so  that 


Eiki<, 

is  proportional  to  n2log  n  log 

Proof:  By  Theorem  5,  the  i’th  residuals  generated  by  both  preconditioned  methods  satisfy 


(24) 


i-i 


<fV, 


.-I 


IMs 

H-oIIq 

where  7<1  is  is  a  constant  that  is  independent  of  mesh  size,  and  q=2  for  CGN,  i)=*l  for 
Orthomin(k).  Hence  (ignoring  the  low  order  contribution  of  *7*=2),  i  >  log(<'l)/log(l/7) 
iterations  suffice  to  reduce  the  relative  error  of  (24)  to  at  most  t.  Since  the  dominant  cost  of 
each  iteration  is  the  preconditioning  solve,  the  total  operation  count  is  proportional  to 
n2log  n  log  <'*. 


Q.E.D. 


The  previous  two  results  are  expressed  in  terms  of  the  Q  -norm  of  the  residuals.  It  might  be  more 
practical  to  measure  the  relative  Euclidean  norm  ||r.||2/||r0||2.  Alternatively,  one  might  wish  to  reduce 
the  error  Sj  ==  u-Uj  to  within  truncation  error.  Since  the  truncation  error  of  the  discretization  is  0(n"2),  it 
is  sufficient  to  reduce  the  relative  error  ||3;||2/||s0||2  by  a  factor  of  0(n'2)  to  compute  a  solution  u;  that  is 
accurate  to  within  truncation  error. 

Corollary  7:  The  number  of  arithmetic  operations  required  by  either  CGN  or  Orthomin(k) 
with  preconditioning  by  Q  to  solve  (I)  so  that 


is  proportional  to  (n2log  n)(log  n  +  log  t*1).  The  number  of  arithmetic  operations  required  to 
compute  a  solution  accurate  to  within  truncation  error  is  proportional  to  n2(log  n)2. 

Proof:  For  the  first  assertion,  note  that 


llroll2  IMq-. 


Since  fc(Q)  =  0(n2)  [10],  the  left  side  of  (25)  is  bounded  by  e  if 


llr0llQ-» 


s<,-  °  U 


The  result  then  follows  from  Corollary  6. 

For  the  second  assertion,  note  that  the  errors  and  residuals  are  related  by  s; 
Hence, 


11  11  —  '  '  11  ti  —  '  '  w  ti  11 


To  reduce  the  relative  error  by  a  factor  of  n’2,  it  suffices  to  make 


*W*(Q) 


Bat,  by  (8)  and  (9), 

<c(A)  <  #e(M)  +  HR)  . 

KJW 


With  M  defined  as  in  (12),  Xmin(M)  —  0(n*2)  and  Xmwt(M)  —  0(1),  so  that  ic(M)  «*  0(n2)  (10). 
Moreover,  as  defined  in  (13),  ail  the  nonzero  entries  of  R  lie  in  four  off-diagonal  bands,  and 
their  absolute  values  are  bounded  by  Ch,  where 

C  =  max  [  max  |e(x,y)|,  max  |d(x,y)|  1 
l  (*j)en  J 

is  independent  of  h  =»  l/(n+l).  By  Gerschgorin’s  theorem  [20],  p(R)  —  0(n'').  Therefore, 
/>(R)/Xmifi(M)  ==  O(n),  so  that  <t(A)  *=*  0(n2).  Substituting  these  results  into  (28), 

<2  0(n'5),  log  tj  =  0(log  »). 

The  result  again  follows  from  Corollary  8. 

Q.E.D. 

The  symmetric  formulation  of  the  preconditioned  problem  (18)  requires  that  the  preconditioning 
matrix  Q  be  symmetric  positive-definite.  For  non-self-ad  joint  problems,  it  seems  preferable  to  allow  Q 
to  be  nonsymmetric  by  including  in  it  a  discrete  separable  approximation  to  the  first  order  terms  in  (2). 
Such  preconditioning  matrices  can  be  used  in  either  of  the  alternative  preconditioned  problems 

AQ^u  =  f,  u  =*  Q^u,  (27) 

Q*‘Au  —  Q'*f,  (28) 

and  the  cyclic  reduction  method  can  be  used  for  the  preconditioning  solves  [17].  An  additional 
advantage  of  (27)  is  that  the  residual  norm  minimized  by  the  two  iterative  methods  under  consideration 
is 

||f  -  AQ-'u^j  =  ||f- Au-|2  , 

i.e.,  it  is  independent  of  the  preconditioning.  However,  we  have  been  unable  to  extend  the  convergence 
analysis  to  these  alternative  problems  (even  for  symmetric  Q).  In  the  analysis  above,  the  application  of 
Theorems  1  and  2  requires  bounds  on  the  extreme  eigenvalues  of  the  symmetric  and  skew-symmetric 
parts  of  A  in  (18).  For  (27),  the  symmetric  part  is  given  by  (AQ'1+[AQ',]T)/2  and  the  skew-symmetric 
part  by  (AQ'1-[AQ',]T)/2.  Wc  have  not  been  able  to  bound  the  eigenvalues  of  these  matrices  or  those 
corresponding  to  (28)  in  a  manner  analogous  to  (18)  and  (23). 


For  an  alternative  approach  to  nonseparable  M,  see  [11]. 


4.  Numerical  Experiments 

In  this  section,  we  present  numerical  results  that  confirm  the  convergence  analysis  of  Section  3  for  (16) 
and  suggest  that  similar  behavior  is  exhibited  for  (27).  All  tests  were  run  on  a  VAX1 1-730  in  double 
precision  (55  bit  mantissa).  The  fast  direct  preconditioning  solves  were  performed  using  the  cyclic 
reduction  method  implemented  in  the  routine  BLKTRI  in  the  FISHPACK  subroutine  package  [19].  We 
consider  two  Dirichlet  problems  of  the  form  (2)  and  one  problem  with  mixed  boundary  conditions. 

For  the  first  two  problems,  let  the  coefficients  of  (3)  be  given  by 

a(x,y)  =  e'xy,  b(x,y)  =»  exy,  c(x,y)  —  0,  (29) 

d(x,y)  =*  l(x+y),  e(x,y)  —  , 

where  7  is  a  scalar  parameter.  The  operator  A  is  nonseparable  and,  for  77^0,  nonsym metric.  The  right 
hand  side  is  determined  by  choosing  the  solution 

u(x,y)  =  x  exy  sin  (jtx)  sin  (sy). 

We  pose  the  problem  on  the  unit  square  0<x,y<l  with  homogeneous  Dirichlet  boundary  conditions,  and 
we  discretize  using  the  five-point  second  order  centered  finite  difference  scheme  on  a  uniform  nin  grid, 
with  h=l/(n+l).  We  use  the  values  7=5,  7=50  and  h  =  1/16,  1/32,  1/64,  1/128,  and  for  one  test, 
1/256. 

We  consider  both  self-adjoint  and  non-self-adjoint  separable  approximations  Q,  which  give  rise  to 
symmetric  and  nonsymmetric  preconditioning  matrices,  respectively.  For  the  self-adjoint  approximation, 
the  coefficients  of  (29)  are  approximated  by 

a(x,y)  —  a(x,5),  b(x,y)  =  b(.5,y),  c(x,y)  =  0, 

d(x,y)  —  0,  e(x,y)  =  ^e(x,.5)+le(.5,y). 

For  the  non-self-adjoint  approximation,  d(y)=d(.5,y).  We  examine  formulation  (16)  with  symmetric  Q 
(the  only  possible  choice),  and  formulation  (27)  with  both  symmetric  and  nonsymmetric  Q.  The 
stopping  criterion  in  all  tests  is 
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HrillQ-*  for  (1®)'  llrill2  f°r(27)- 

The  initial  guess  is  u0=0. 


Tables  4*1  and  4*2  show  the  iteration  counts  for  7=5  and  7=50,  respectively.  In  the  tables,  both  the 
preconditioning  formulation  and  its  associated  norm  are  listed. 


To  examine  different  boundary  conditions,  for  the  third  problem  we  consider  the  differential  equation 

[12] 


-(U„  +  «yy)  +  “  0 

on  {(x,y)  |  x>0,  y>0),  with  boundary  conditions 


u(x,0)  =  0,  u(0,y)  =  1, 
u(x,y)  bounded  as  |x|  +  |y|  -*  oo. 


The  exact  solution  to  (30)  -  (31)  has  a  boundary  layer  at  y=0  and  is  nearly  identically  one  elsewhere. 
Following  [12],  for  the  discrete  problem  we  restrict  (30)  to  the  unit  square  and  add  the  boundary 


conditions 


u(x,l)  =  l,  ux(l,y)  =  0.  (32) 

We  discretize  on  a  uniform  nin  grid  using  centered  differences  for  all  terms  except  at  the  right  boundary, 
where  we  use  the  first  order  approximation 


U  ,  .  •  U  • 
_n+lj__ni 


0  ~  ~  \  (w) 

Incorporating  (33)  into  the  equations  centered  at  {u  •}?  .  results  in  a  linear  system  of  order  N=n2.  We 

uj  j  ■ 


consider  the  value  /?=10. 


Because  (30)-(32)  is  separable,  the  discrete  problem  can  be  solved  directly  by  cyclic  reduction.  We 
therefore  consider  only  symmetric  preconditioning  matrices  based  on  the  discrete  Laplacian,  with  the 
discretization  at  the  right  boundary  handled  as  in  (33).  The  iteration  counts  for  the  same  stopping 
criterion,  mesh  sizes  and  initial  guess  as  above  are  shown  in  Table  (4-3). 
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Table  4-1:  Discretization  of  equation  (29),  7=5.  Iterations  to  reduce  residual 

norm  by  factor  of  10‘®. 
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Table  4-2:  Discretization  of  equation  (29),  7=50.  Iterations  to  reduce  residual 

norm  by  factor  of  10'®. 


*For  a  given  non-self-adjoint  operator,  eydie  reduction  is  applicable  only  for  small  enough  h  (see  | I7J,  p.  IH2,  or  1 10)).  For 
7—50,  h  —  1/10  and  1/32  are  too  large.  These  mesh  sixes  are  handled  in  the  experiments  by  reducing  the  contribution  to  Q 
of  the  first  derivatives,  taking  <T(y)— i(y)d(.5jr)  with  0<6(y)<  1,  so  the  approximation  of  R  in  Q  is  probably  less  accurate  in 

these  eases. 
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Table  4-3:  Discretization  of  equation  (30)-(32)>  ^=10.  Iterations  to  reduce  residual 

norm  by  factor  of  10‘®. 

These  experiments  confirm  the  convergence  analysis  of  (16)  and,  with  one  exception,  suggest  that 
convergence  for  (27)  is  independent  of  mesh  size  also,  for  both  symmetric  and  nonsymmetric  Q.  The 
exception  occurs  with  7=50,  for  which  the  iteration  count  for  CGN  is  still  growing  at  h=  1/256,  and 
Orthomin(l)  fails  to  converge  for  h<  1/128,  indicating  that  (27)  may  require  a  very  fine  mesh  if  A  is  not 
well-approximated  by  a  symmetric  preconditioning  matrix.  (The  convergence  failures  are  due  to  the  fact 
that  the  symmetric  part  (AQ',+(AQ'1)T)/2  of  the  coefficient  matrix  of  (27)  is  indefinite  for  the  four 
mesh  sizes  considered.)  The  smaller  iteration  counts  in  Tables  4-1  and  4-2  for  nonsymmetric 
preconditioning  matrices  suggest  that  nonsymmetric  preconditioners  offer  some  advantage.  However,  the 
operation  counts  for  cyclic  reduction  on  nonsymmetric  matrices  (20n2log  n  +  0(n2))  are  higher  than 
those  for  fast  direct  methods  for  symmetric  matrices  (e.g.,  5n2log  n  +  0(n2)  for  cyclic  reduction;  see 
(2,  3,  7,  18]).  For  the  two  problems  examined,  the  asymptotic  counts  for  symmetrically  applied 
symmetric  preconditioners  are  actually  lower. 

Finally,  we  note  that  a  comparison  between  these  preconditioning  techniques  and  others  based  on 
incomplete  factorizations  can  be  found  in  [9]. 
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